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ABSTRACT 



Dust particles sediment toward the midplanes of protoplanetary disks, forming dust-rich sub- 
layers encased in gas. What densities must the particle sublayer attain before it can fragment 
by self-gravity? We describe various candidate threshold densities. One of these is the Roche 
density, which is that required for a strengthless satellite to resist tidal disruption by its pri- 
mary. Another is the Toomre density, which is that required for de-stabilizing self-gravity to 
defeat the stabilizing influences of pressure and rotation. We show that for sublayers containing 
aerodynamically well-coupled dust, the Toomre density exceeds the Roche density by many (up 
to about 4) orders of magnitude. We present 3D shearing box simulations of self-gravitating, 
stratified, dust-gas mixtures to test which of the candidate thresholds is relevant for collapse. 
All our simulations indicate that the larger Toomre density is required for collapse. This result 
is sensible because sublayers are readily stabilized by pressure. Sound-crossing times for thin 
layers are easily shorter than free-fall times, and the effective sound speed in dust-gas suspen- 
sions decreases only weakly with the dust-to-gas ratio (as the inverse square root). Our findings 
assume that particles are small enough that their stopping times in gas are shorter than all other 
timescales. Relaxing this assumption may lower the threshold for gravitational collapse back 
down to the Roche criterion. In particular, if the particle stopping time becomes longer than the 
sound-crossing time, sublayers may lose pressure support and become gravitationally unstable. 

Subject headings: hydrodynamics — instabilities — planets and satellites: formation — protoplanetary 
disks — methods: numerical 



Gravitational instability is an attractive mech- 
anism to form planetesimals, but how it is trig- 
gered in protoplanetary disks remains unclear. In 
one proposed sequence of events, most of the disk's 
solids first coagulate into particles 0.1-1 m in size 
at orbital distances of a few AU. These "boulder"- 
sized bodies then further concentrate by the aero- 
dynamic streaming instability (Youdin & Good- 
man 2005; Johansen et al. 2007; Bai & Stone 
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2010; and references therein). Local densities are 
so strongly enhanced by the streaming instability 
that they can exceed the Roche density (see §1.3 
for a definition), whereupon collections of boul- 
ders may undergo gravitational collapse into more 
massive, bound structures. 



A weakness of this scenario is that it presumes 
that particle-particle sticking (i.e., chemical ad- 
hesion) can convert most of the disk's solids into 
boulders, or more accurately, particles whose mo- 
mentum stopping times in gas 
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are within a factor of 10 of the local dynamical 
time f2 _1 , where fi is the Kepler orbital frequency, 
to is the particle mass, v rc \ is the relative gas- 
particle velocity, and -fdrag is the drag force whose 
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form varies with disk environment (see, e.g., Wei- 
dcnschilling 1977). Figure 1 relates i st0 p to parti- 
cle radius s as a function of disk radius r in the 
minimum-mass solar nebula (MMSN). For r ps 1- 
10 AU, the condition £lt stop — 0.1-1 corresponds 
to s ps 0.1-1 m. 




Disk radius r (AU) 

Fig. 1. — Stopping times of particles in the MMSN, 
normalized to the local dynamical time Q~ . Disk 
parameters are taken from Chiang & Youdin (2010). 
Particles are assumed spherical with bulk density 1 
g/cm 3 . The kinks in the curves are due to transitions 
between different drag force laws as taken from Wei- 
denschilling (1977; note that the transition between 
the Stokes and Epstein drag laws occurs when the gas 
mean free path equals 2/9 of the particle radius, not 
4/9 as misprinted in that article). Marginally coupled 
particles (fit s t op ~ 1) correspond to meter-sized boul- 
ders at r ~ 1 AU; decimeter-sized rocks at r ~ 10 
AU; and cm-sized pebbles at r ~ 100 AU. The top 
panel plots f2t s to P at various fixed particle radii s; the 
bottom panel plots the same data but at fixed Qt a t op - 
In this paper we are interested in the small particle 
nt s top <C 1 limit. 

Unfortunately, particle-particle sticking might 
not produce boulders in sufficient numbers for 
the streaming instability to be significant. A 



comprehensive study by Zsom et al. (2011; see 
also Birnsticl ct al. 2010) found that for real- 
istic, experiment-based sticking models that in- 
clude both bouncing and fragmentation, particles 
no larger than ~1 cm can form by sticking — even 
when the disk is assumed to have zero turbulence. 
According to Table 1 of Zsom et al. (2011), coag- 
ulation models over most of parameter space pro- 
duce t s ~ 10 -4 -10 -2 . This range is too small 
for the streaming instability to concentrate parti- 
cles strongly — see Bai & Stone (2010), who showed 
that when half or more of the disk's solid mass has 
f^stop < 0.1, densities enhanced by the streaming 
instability still fall short of the Roche density by 
more than a factor of 10. Even if particle-particle 
sticking could grow bodies with £lt s t op ~ 0.1-1 
(e.g., Okuzumi et al. 2009, who neglected fragmen- 
tation), the disk's solids may not be transformed 
into such bodies all at once. Rather, boulders may 
initially comprise a minority on the extreme tail 
of the size distribution. Unless they can multiply 
from a minority to a majority within the time it 
takes for them to drift radially inward by gas drag 
(-100-1000 yr starting at 1 AU; Wcidenschilling 
1977), they threaten to be lost from the nebula by 
drag. 

We are therefore motivated to ask whether 
gravitational instability is practicable for parti- 
cles having realistically small sizes and concomi- 
tantly short stopping times, say Oi stop < 10 -2 . 
Smaller particles suffer the disadvantage that they 
are harder to concentrate; since they are well- 
entrained in gas, turbulence in the gas can loft 
particles above the midplane and prevent them 
from collecting into regions of higher density. The 
streaming instability provides one source of turbu- 
lence. Another driver of turbulence is the Kelvin- 
Helmholtz instability, caused by vertical velocity 
gradients which steepen as dust settles into a thin, 
dense "sublayer" at the disk midplane (Weiden- 
schilling 1980). Several recent studies (Chiang 
2008, Lee et al. 2010a,b; see also Weidenschilling 
2006, 2010) have measured the maximum sublayer 
densities permitted by the Kelvin-Hclmholtz in- 
stability. Neglecting self-gravity, they found that 
dust-to-gas ratios between ~2-30 are possible in 
disks that are locally enriched in metallicity by fac- 
tors of 1-4 above solar. Such local enrichment can 
be generated by radial drifts of particles relative 
to gas (see Chiang & Youdin 2010 for a review). 
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For observational evidence of radial segregation of 
dust from gas, see Andrews et al. (2012). 

Are such enhancements in the local dust-to-gas 
ratio sufficient to spawn planetesimals? How high 
must dust + gas densities be before the effects of 
self-gravity manifest? Our paper addresses these 
questions in the limit Qt stop <C 1, i.e., in the 
limit that particles are small enough to be well 
coupled to gas. In the next two subsections, we 
derive critical densities for gravitational instabil- 
ity in the cases of a pure gas disk (§1.1), and a 
disk composed of both gas and perfectly entrained 
(fiistop — > 0) dust (§1-2). The two cases give re- 
markably different answers for dust-rich sublayers. 
In §1.3 we add two more densities from the liter- 
ature to the list of proposed criteria for gravita- 
tional collapse. Table f summarizes the various 
candidate threshold densities. 

In §2-§3, we present numerical simulations of 
3D, self-gravitating, compressible flows of thin, 
dense sublayers of dust. We use these simulations 
to try to identify which of the proposed criteria (if 
any) is the most relevant for gravitational insta- 
bility. Section 4 summarizes our findings but also 
points out the limitations of our numerical simula- 
tions, which are restricted to the asymptotic limit 
f^stop — > 0. We argue in §4.1 how finite but still 
small values of H,t stop may lower the threshold for 
gravitational collapse. 

1.1. Critical Density for Gravitational 
Instability in a Pure Gas Disk 

The usual criterion for gravitational instability 
in a razor-thin pure gas disk is expressed in terms 
of the dimensionless parameter 



Q g 



ttGXL 



(2) 



where G is the gravitational constant, c g is the 
gas sound speed, and £ g is the gas surface den- 
sity (Goldreich & Lynden-Bell 1965; Toomre 1964; 
Toomrc 1981). In (2), the Kepler orbital frequency 
£1 has been substituted for the radial epicyclic fre- 
quency. If 

Q g < Qg - i , (3) 

the disk is gravitationally unstable to axisym- 
metric perturbations in the disk plane. The 
Q-criterion is a measure of the competition be- 
tween stabilizing pressure, stabilizing rotation, 



and de-stabilizing self-gravity (see, e.g., Binney 
& Trcmainc 2008). When Q g > 1, horizontal 
perturbations having lengthscales < 2c g /G£ g 
are stabilized by pressure, while those having 
lengthscales > 2c g /GE g are stabilized by rota- 
tion. When Q g equals 1, the first axisymmctric 
mode to become unstable to self-gravity has ra- 
dial wavelength 2c g /G£ g . And as Q g approaches 
1 from above, the disk is increasingly susceptible 
to nonaxisymmetric perturbations which swing 
amplify (Goldreich & Lyndcn-Bcll 1965). 

The criterion Q g < Q* for gravitational insta- 
bility can be translated into a criterion for the 
midplane density p g o (the subscript "0" denotes 
the initial midplane value) . We define a disk half- 
thickness H g using 

E g = 2pg H g . (4) 

We also define a half-thickness using the usual 
relation from vertical hydrostatic equilibrium: 



cg/n. 



(5) 



Ordinarily H g ss ifj and we would not bother to 
distinguish the two; however, we will find cases 
later on where they differ by factors of several be- 
cause of the effects of dust, and thus we take care 
to separate the two lengths now. Upon substitu- 
tion of (4) and (5), the relation Q g < Q* is shown 
to be equivalent to 

i m + 



Pgo > P*i 



1 



P' 



(6) 



2tt Q* Hg 

where we have defined a reference density 1 

P f = Mjr 3 (7) 

with and r equal to the mass of the central 
star and the disk radius, respectively. 

The pj-criterion (6) is sometimes used (e.g., Lee 
et al. 2010a, b) to signal gravitational instability in 
dusty gas disks (with p g0 replaced by the total dust 
-I- gas density p da + p g0 , Q* R = 1, and fl|/-Hg = 1). 
But using p\ for dust-gas mixtures is suspect be- 
cause the criterion does not account explicitly for 
the two-phase nature of such media. In the next 
subsection we make such an accounting to derive 
a substantially different criterion for gravitational 
collapse. 



1 In this paper, we will superscript critical threshold densities 
with *, and fiducial or reference quantities with f. 
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Table 1 

Candidate Critical Densities for Gravitational Collapse 



Critical 
Density 




Value 


Comment 


Reference 


pi 


l i 


^i% -0.16% W 

rig r° r° 


Equivalent to Q g < Q* for pure gas 
disks 


This paper, equation (6) 


PSekiya 




0.60% 


Required for the onset of an in- 
compressible, axisymmctric ovcrstablc 
mode 


Sekiya (1983) 


''Roche 




3.5% 


Required by satellite to resist tidal dis- 
ruption by primary 


Chandrasekhar (1987) 




1 Qg /Eg 


) 2 g%-2xl0 4 % 


Equivalent to < Q*j for dust-rich 
sublayers in gas 


This paper, equation (13) 


(«) Value 


is derived for Q 


'* = 1 and Hl/H s = 1. 






( 6 ) Value 


is derived for Q 


* = 1, Q g = 30, S d /S g 


= 0.015, and H^/Hg = 1. 





1.2. 



Critical Density for Gravitational 
Instability in a Dust-Rich Sublayer 



in the Limit fit 



stop 







For disks of gas and dust, gravitational instabil- 
ity should still be determined by the Q-criterion, 
except there is now the possibility that disk self- 
gravity is dominated by dust in a vertically thin 
sublayer at the midplane: 



Qc, 



< 



Q*j for instability. 



(8) 



In using the dust surface density Ed in (8), we 
neglect the contribution of gas to the total sur- 
face density of the sublayer. Under typical cir- 
cumstances, the error accrued is small. 



In the limit fit 



stop 



0, the dust-gas mixture 
represents a colloidal suspension. In this suspen- 
sion, dust does not contribute to the pressure P 
— which is still provided entirely by gas — but 
instead adds to the inertia. In other words, 



P : 



Pd)c\ 



(9) 



by definition of Cd, the speed of sound in the sus- 
pension: 



V 1 + P 

where p g is the local gas density, pd is the local 
dust density, and p = Pd/Pg is the dust-to-gas 
ratio. In effect, dust increases the mean molecular 
weight of the gas. 



Inserting (10) into (8) and using 

= t 1 H t 
s 2nQg H g 



(11) 



we solve for the total midplane density required 
for gravitational instability: 



where 



Pq = Pdo + Pgo <; Pii 



Pn 2ir Q* 2 V Ed J H. 



(12) 



(13) 



4 A I Qg\ ( 1 



2 x l(Tp , 

30 ) \Q* d 
0.015 \ 2 / Hl/Hg 



(14) 



In (14), our normalizations for Q g and the bulk 
(height-integrated but local to r) metallicity 
E d /E g derive from the MMSN at r = 1 AU (Chi- 
ang & Youdin 2010). For these parameter choices, 
the critical midplane density pjj is an astonish- 
ing five orders of magnitude greater than p\. It 
is possible that real disks have masses and bulk 
metallicities enhanced over the MMSN by factors 
of a few, in which case pj : would be larger than p\ 
by about three orders of magnitude. 
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1.3. Other Critical Densities 



2. METHODS 



Another threshold density, already alluded to 
at the beginning of §1, is the Roche density: 

PRoche = 3-5— jj- • (15) 

The Roche density is the density required for a 
strengthless, incompressible, fluid body in hydro- 
static equilibrium to resist tidal disruption, when 
in synchronous orbit at distance r about a star of 
mass (e.g., Chandrasekhar 1987). 

Yet another candidate threshold was proposed 
by Sekiya (1983), who found that when the mid- 
plane density exceeds 

PSckiya = °- 60 ^f » ( 16 ) 

the disk becomes susceptible to an unstable, in- 
compressible, axisymmetric mode in which in- 
plane motions generate out-of-plane bulges (i.e., 
an annulus that contracts radially becomes thicker 
vertically, and vice versa). The nonlinear outcome 
of this instability is not known, but Sekiya (1983) 
speculated that the dust sublayer might eventu- 
ally fragment on the scale of the wavelength of 
the overstable mode, and that dust particles might 
sediment toward the centers of fragments to form 
the first-generation planetesimals. 

Table 1 summarizes the four candidate thresh- 
old densities. For realistic parameters (Q g ~ 10- 
30; £ d /£ g ~ 0.015-0.15), the four densities obey 

Pi < PSckiya < PRoche < Pu ■ ( 17 ) 

The smallest three densities in this hierarchy 
are fixed multiples of the reference density = 
M«/r 3 (with coefficients ~1/2tt w 0.16, 0.6, and 
3.5, respectively). The last density pjj can, in 
principle, be arbitrarily larger than p'; for typi- 
cal, astrophysically plausible parameters, it is 2-4 
orders of magnitude larger. 

Which of the four densities in Table 1 is the 
most accurate predictor of gravitational collapse? 
In the next two sections, we describe numerical 
simulations performed in the ftt s top — > limit that 
attempt to answer this question. We will find un- 
fortunately that the numerical expense of simulat- 
ing thin sublayers of dusty gas will force us into a 
parameter space where the difference between pjj 
and the other densities is not as large as it is in 
reality; we will have to make do with what we can. 



2.1. Code 

We simulate hydrodynamic, self-gravitating, 
stratified flows in disks using Athena, configured 
for a shearing box, with no magnetic fields (Stone 
et al. 2008; Stone & Gardiner 2010). Dust is as- 
sumed to be perfectly aerodynamically coupled to 
gas so that they share the same velocity field v. 

The equations solved are: 

g+V-(pv) = 0, (18) 
^+V-( Pg v) = 0, (19) 

^ + V • (pvv + P) = -pV$ 

-2/?(fiz) x v + 2qpD, 2 x± - ptt 2 zz , (20) 
V 2 $ = 4ttGp, (21) 

where p — p g + is the total density of the dust- 
gas suspension, P = PI is a diagonal tensor with 
components P — PgC 2 as defined in equation (9) 
with constant c g (isothermal approximation), Q is 
the mean (constant) orbital frequency, x points in 
the radial direction, z points in the vertical direc- 
tion, and <& is the self-gravitational potential of the 
dust-gas mixture. We choose the shear parameter 
q = 3/2 for Keplerian flow. 

2.1.1. Algorithms and boundary conditions 

Athena 4 . provides several schemes for time 
integration and spatial reconstruction, and for 
solving the Riemann problem. Having experi- 
mented with various options, we adopted the van 
Leer algorithm for our dimensionally unsplit inte- 
grator (van Leer 2006; Stone & Gardiner 2009); a 
piecewise linear spatial reconstruction in the prim- 
itive variables; and the HLLC (Harten-Lax-van 
Leer-Contact) Riemann solver. To account for 
disk self-gravity, we use the routines written by 
Koyama & Ostriker (2009) and Kim et al. (2011) 
which solve Poisson's equation using fast Fourier 
transforms. 

Boundary conditions for our hydrodynamic 
flow variables (including density and velocity, but 
not the self-gravitational potential) are shearing- 
periodic in radius (x) and periodic in azimuth 
(y). For vertical height (z), we experimented with 
both periodic and outflow boundary conditions, 
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and chose periodic boundary conditions to ensure 
strict mass conservation. When outflow boundary 
conditions were employed, mass was lost from the 
boundaries at early times and complicated the in- 
terpretation of our results. We verified that our 
results are insensitive to box height for sufficiently 
tall boxes; see §3 for explicit tests. 

The Poisson solver implements shear-periodic 
boundary conditions in x, periodic boundary con- 
ditions in y, and vacuum boundary conditions in 
z (Koyama & Ostriker 2009; Kim et al. 2011). In 
our simulations, self-gravity is dominated by dust, 
and our boxes are tall enough to contain the en- 
tire dust layer. Both vertical and radial stellar 
tidal gravity are included as source terms in the 
van Leer integrator. 

We further augmented the code to include dust 
in the limit of zero stopping time. In this limit, 
dust shares the same velocity field as gas, and con- 
tributes only to the mass density. In our modified 
version of Athena, two continuity equations are 
solved: one for the entire mixture (p = p g + pd, 
see equation 18), and one for the gas (p g , see equa- 
tion 19). The dust density is given by the dif- 
ference (p — p g ). The remaining hydrodynamic 
equations govern the dust-gas mixture (p), but 
with gas (p g ) contributing solely to the pressure P 
(see Equation (20)). For simplicity, we adopt an 
isothermal equation of state so that P is related 
to p g by equation (9) for constant c g . Isothermal 
flows are more prone to gravitational instability 
than adiabatic ones (Mamatsashvili & Rice 2010). 

We also modified the HLLC Riemann solver to 
accommodate our dust-gas mixture. Changes in- 
clude the following: (1) The speeds of the left, 
right, and contact waves are reduced by a factor 
(1 + p)~ x / 2 , where p = Pd/p g is the local dust- 
to-gas ratio, to account for the added inertia from 
dust (see equation 10). (2) The pressure in the 
contact region is replaced by an equivalent but 
numerically more accurate form based on equa- 
tion (10.76) in Toro (1999). (3) When calculating 
left/right momentum fluxes, we ensure that only 
gas contributes to the pressure by using p g and 
not p. (4) For the flux solver to predict the pres- 
sure and wave speeds, the left /right gas densities 
require specification. We therefore add a recon- 
struction process for the gas density which inter- 
polates cell-centered values to cell boundaries to 
second-order accuracy. 



Previous studies of dust in the perfectly coupled 
limit (Chiang 2008; Lee et al. 2010a, b) also intro- 
duced a static background radial pressure gradi- 
ent to mimic sub-Keplerian rotation of gas in a 
pressure-supported disk. We could also add the 
appropriate source term to the van Leer integra- 
tor. However, since our goal is to determine the 
minimum densities required for gravitational col- 
lapse and not to study vertical shear instabilities 
(i.e., the Kelvin-Helmholtz instability), we omit 
the background pressure gradient for simplicity. 

In many of our simulations, the dust layer at the 
midplane collapses vertically because it is gravita- 
tionally unstable. Because of our boundary condi- 
tions, "fresh" gas from outside the simulation box 
cannot enter into the box, and thus in the event 
of gravitational collapse toward the midplane, the 
topmost and bottommost regions of our simula- 
tion domain become evacuated. Low-density gas 
in those regions become increasingly easy to ac- 
celerate, and the code timestep shortens by orders 
of magnitude, effectively halting the simulation. 
The dramatic reduction in timestep is not a seri- 
ous limitation, as it usually occurs after the col- 
lapsing dust has attained some saturated state (see 
§3.2.1). In any case, we are more interested in the 
onset of gravitational instability than its nonlinear 
development. 

2.1.2. Code tests 

The following test problems helped to validate 
our code. 

Linear wave propagation. — We propagated a 
small-amplitude ID wave in a medium with a uni- 
form background dust-to-gas ratio, with periodic 
boundary conditions, no background shear, and 
no gravity. We found the simulated wave speed 
matched the reduced sound speed calculated in 
(10). We chose our box to be one wavelength long, 
so that after one wave period, the wave crossed 
the boundaries and returned to its original posi- 
tion. With N = 128 grid cells and an initial (frac- 
tional) wave amplitude A = 10~ 4 , we found the 

N 

deviation Sq = j^/Jl^i — qi\ ~ 2 x 10~ 8 , where 

i=l 

q G {pd, Pg, p} and q° represents the initial condi- 
tion. 

Dust cloud advection. — We advected a Gaussian- 
shaped dust cloud in a ID domain. The cloud oc- 



G 



cupied about half the size of the box and the code 
was run for one box-crossing time. With N = 256 
grid cells, the root-mean-squared deviation in the 
shape of the cloud was < 1%. 

Hydrostatic equilibrium of a stratified but non- 
self-gravitating dusty disk. — Omitting self-gravity 
but including stellar gravity (both radial and ver- 
tical), we set up 3D dust-gas mixtures in hydro- 
static equilibrium. A variety of vertical profiles 
for the dust-to-gas ratio were tested, ranging from 
uniform to linear to more complicated functional 
forms. All equilibria were found to be stable 
against small perturbations, even for dust-to-gas 
ratios as large as several hundred. 

Gravitational instability of 3D pure gas disks. — 
We simulated isothermal, gravitationally unstable 
disks of pure gas in 3D. The gas was initialized 
in hydrostatic equilibrium (computed with verti- 
cal self-gravity), and box heights spanned approx- 
imately ±4 initial gas scale heights. We found 
that Q g = 1 did not trigger gravitational insta- 
bility, whereas Q g = 0.5 did. Our results are 
consistent with those of Goldrcich & Lyndcn-Bcll 
(1965), who found analytically that Q* = 0.676 
for a finite-thickness isothermal gas disk. 

2.2. Initial Conditions 

and Run Parameters 

Initial conditions for our science simulations 
are of a dust-gas mixture with a pre-defined 
vertical profile for the dust-to-gas ratio p(z) = 
Pa(z)l 'pg(z). We choose the form 



Taking derivatives, we find 



p{z) = p sech ' 



(22) 



where po is the midplane dust-to-gas ratio. The 
scale height za can be thought of as the half- 
thickness of the dust layer insofar as p g (z) is con- 
stant with z. 

The isothermal dust-gas mixture is initialized 
in vertical hydrostatic equilibrium, including both 
stellar tidal gravity and disk self-gravity: 

^_^ = _ fi 2 z _ 47rG [\ +pd)dz . (23) 
p g + pa dz J 

We solve numerically the differential form of (23). 



d 
dz 



(1 + P) 



-ldlnpg 



dz 



Hi 



4vrG 



(24) 

where iJ| = c g /Q is a fiducial (constant) gas scale 
height, not to be confused with any actual disk 
scale height. A non-dimensional form of (24) is 
given by 



d 

dz 



(1 + P) 



dz 



h g Q, 



Pg^ + P), 



. (25) 

where we have defined the dimensionless variables 
z = z/H^, p g = Pg/p g o (where p g o is the midplane 
gas density), and h g = H g /H^, with 



"g — -"g/ 
H g = S g /(2p g0 ) . 



(26) 



Upon insertion of (22), equation (25) can be 
solved numerically for p g (z). But the solution 
must satisfy the following two constraints: 



h s = / Pg(z)dz 



by definition of H g , and 



Ed 
E„ 



1 



,{z)n{z)dz 



(27) 



(28) 



for a fixed height-integrated (i.e., bulk) metallicity 
Ed/S g . 

Our procedure is as follows. We freely specify 
Q g , Ed/Eg, and po as model input parameters. 
We then iteratively solve equations (25), (27) and 
(28) for the three unknowns /5 g (z), h s , and z^. 
First we guess Zd and h g , and integrate (25) to 
obtain p g (z). If p g so calculated fails (27), then 
we revise h g and re- integrate (25), repeating un- 
til (27) is satisfied. Next we check (28). If p s (z) 
and hg fail (28), then we revise Zd and repeat the 
procedure from the beginning, re-integrating (25) 
to obtain p g , re-establishing (27), and so on. Typ- 
ically ^100 iterations (~10 for z& x ~10 for h g ) 
are required before all constraints are satisfied to 
~1% accuracy in z^ and 10 -6 accuracy in hg. 

Table 3 lists the parameters of our models. 
Note that these parameters do not describe plau- 
sible protoplanetary gas disks; in particular, our 
model metallicities Sd/E g are orders of magnitude 
above the solar value of ^0.015. Parameters are 
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instead chosen to yield disk flows that our code 
can adequately resolve while still testing equation 
(14). Unfortunately, more astrophysically realis- 
tic parameters correspond to dust sublayers that 
are too vertically thin for us to resolve numeri- 
cally; the code timestep, set by the sound-crossing 
time across a grid cell, becomes prohibitively short 
as thinner dust layers are considered. This diffi- 
culty means that the difference between pjj and 
the other candidate threshold densities is much 
less than what it is in reality and our ability to 
distinguish between the candidates degrades as a 
result. 

Figure 2 plots the initial conditions for our stan- 
dard model (S = STD32), for which Q g = 24, 
Sd/Sg = 8, and /io = 35. For this specific 
case, we calculate that h g = 0.20 c g /i! and z& = 
0.083 Cg/n. The top and bottom boundaries of 
our simulation box are indicated by dotted ver- 
tical lines; typically box heights span ±4z<i (see 
§3 for box height tests). The right-hand panel of 
Figure 2 compares gas density profiles computed 
with and without self-gravity, and with and with- 
out dust, and shows that both the weight and self- 
gravity of the embedded dust layer force gas into 
a similarly thin layer. 




-0.4 -0.2 0.0 0.2 0.4 1 2 3 

z (c a /£l) z (c g «2) 



Fig. 2. — Left: Initial dust and gas densities for 
standard run S = STD32. Right: Initial gas den- 
sity (dashed) on an expanded scale, together with the 
gas density computed without self-gravity but with 
dust (dash-dot) and without dust but with self-gravity 
(dash-double-dot). Vertical lines in both panels lines 
delimit the top and bottom of our computational box. 

Figure 3 plots the initial force densities within 
the upper half of the simulation box to demon- 
strate how well vertical hydrostatic equilibrium is 
satisfied. The sum of stellar gravity (blue dashed 
curve) and disk self-gravity (red solid curve com- 
puted via the integral in equation 23) should equal 



1.5 



a i.o 



o 0.0 




-0.5 r ................................ i 

0.00 0.05 0.10 0.15 0.20 0.25 0.30 

z ( Cfl /n) 

Fig. 3. — Vertical force balance for our initial con- 
ditions. The force density due to self-gravity is com- 
puted two ways: by direct integration of the density 
profile (red solid), and by using the code's Poisson 
solver (green dashed). The two methods agree. The 
force density due to the pressure gradient (black solid) 
should equal the sum of self-gravity and the static stel- 
lar potential (blue dashed). The horizontal gray dot- 
ted line shows the ratio of pressure to gravity. All force 
densities are shown in their absolute values. 

the pressure gradient (black solid curve). It does, 
as evidenced by the ratio of pressure to gravity 
(gray dotted line) which is practically constant at 
unity. We also overplot the self-gravitational force 
computed by our 3D Poisson solver (green dashed 
curve); the agreement with the exact solution is 
good. 

Every simulation listed in Table 3 is perturbed 
from its initial equilibrium by adding random cell- 
to-cell fluctuations of amplitude ~10 _3 c g to the 
velocity field. The typical duration of a simula- 
tion is ^20 f^ 1 . Our rationales for box size and 
resolution are explained in §3. 

3. RESULTS 

Results for 2D shearing sheets are described in 
§3.1, and those for 3D shearing boxes are in §3.2. 

3.1. 2D Shearing Sheet 

For two-dimensional dusty disks, the criterion 
for gravitational instability reads 

Qd ~ irG(V d + Z s )- (1 + M0 )3/2 <y d,2D- 

(29) 
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Table 2 
2D Simulation Parameters 



Name 


Q s 


Mo 


Qd 


A c (cg/fi) 


L x x L y (c 2 g /Q 2 ) 


Resolution 


GlO) 


Duration (S7 ) 


S2D0 


12 


8 


0.44 


0.93 


10 x 10 


256 X 256 


Y 


5.9 


S2D1 


12 


8 


0.44 


0.93 


10 x 10 


512 x 512 


Y 


6 


S2D2 


12 


8 


0.44 


0.93 


0.5 x 0.5 


32 x 32 


N 


100 


S2D3 


12 


8 


0.44 


0.93 


0.5 x 0.5 


128 x 128 


N 


100 


S2D4 


12 


1.2 


1.0 


2.8 


30 x 30 


256 x 256 


N 


100 


S2D5 


12 


2.3 


2.0 


6.9 


70 x 70 


256 x 256 


N 


100 


S2D6 


6 


4.2 


0.5 


1.4 


15 x 15 


256 x 256 


Y 


6.2 



W GI = Gravitational Instability. Y means maxSj increases by orders of magnitude over a few 
dynamical times, and N means it does not. 



We test this criterion by constructing a series of 
2D shearing sheet simulations with various values 
of Q g and Mo , thereby seeing if we can converge on 
a unique value for Q*j 2D . Although total surface 
densities can change during the simulation, the 
dust-to-gas ratio stays fixed at its initial value be- 
cause of our perfect-coupling approximation. Ini- 
tial conditions are as follows: for a given domain 
size L x and L yi the flow velocity v = — ^£lxe y and 
the surface density £ = So + <S£cos(k ■ x), with 
S = S g0 + S d o = Sgo(l + Mo), <5£/£ = °- 01 > 
k x = — 2(2ir/L x ), and k y — 2ir/L y . In our 2D 
simulations, we choose c g = fi = S g o = 1 as our 
units. 

Table 2 lists the parameters for our 2D runs. 
Our standard 2D run, labeled S2D0, has Q g = 12 
and fio = 8.0 and therefore Qd = 0.44. For this 
run, the domain size is chosen large enough to eas- 
ily fit the critical wavelength A c for gravitational 
instability: L x = L y = 10c g /f2 > 10A C , where 

is the wavelength of the fastest growing mode ac- 
cording to the WKB dispersion relation for ax- 
isymmetric waves. It is also the wavelength of 
the first mode to become unstable when Qd just 
crosses Q* A 2D . The resolution of the standard run 
is N x x N y — 256 x 256 so that one critical wave- 
length is resolved across ~10 grid cells. 

For S2D0, we find that the disk is indeed gravi- 
tationally unstable: density waves steepen quickly, 
and dense clumps of dusty gas form before one 
orbital period elapses. A simple way to portray 
instability is to track the maximum dust density 



maxSd versus time — this is done in Figure 4, 
which shows that the maximum dust density in- 
creases by two orders of magnitude over a few dy- 
namical times for our standard run. 



1000 




2 4 6 8 

t (tr 1 ) 



Fig. 4. — Time evolution of the maximum dust den- 
sity max Ed for our 2D shearing sheet simulations. 
The critical value Qd,2D below which gravitational in- 
stability is triggered appears to be between 0.5 and 
1.0. 

Also shown in Figure 4 are results for other 
runs. In S2D4, S2D5, and S2D6, either Q g or 
Mo is varied relative to our standard run, so that 
Qd varies from 0.5 to 2.0. For all these runs, the 
domain size is ^10A C in each direction and the res- 
olution is ~ 10 cells per A c , just as in the standard 
case. Taken together, the results indicate that 

0.5<Q d , 2D < 1.0. (31) 

Other runs explore the effects of varying resolu- 
tion and domain size. Doubling both N x and N y 
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relative to our standard run (S2D1) enables higher 
maximum densities to be achieved when the insta- 
bility saturates, but otherwise does not seem to 
alter the evolution. Reducing the size of the box 
so that it can no longer accommodate even a sin- 
gle critical wavelength (S2D2, S2D3) results in no 
instability, as expected (Gammie 2001; Johnson & 
Gammie 2003). 

3.2. 3D Stratified Dusty Disks 

Equation (8; equivalently 29) gives the crite- 
rion for gravitational instability in a 2D razor-thin 
sheet. For a 3D, vertically stratified disk, there is 
some ambiguity as to how we evaluate Cd in equa- 
tion (8) because its value varies with height. Here 
we simply take Cd to be its value at the midplane, 
so that criterion (8) becomes 

* Q g n * W2 < Q* d ■ (32) 

An alternative is to calculate a vertically averaged, 
density-weighted sound speed. We found, how- 
ever, that such a procedure made little practical 
difference, since dust densities are much greater 
than gas densities near the midplane and drop 
steeply with height. 

3.2.1. Standard run (S = STD32) 

To orient the reader, we present results for our 
standard 3D run (S, also labeled STD32 in §3.2.2), 
for which Qa = 0.5. The full set of model S pa- 
rameters are listed in Table 3, and the initial gas 
and dust density profiles are displayed in Figure 2. 
Our simulation box extends ±4zd vertically, and 
14zd in either horizontal direction. Each horizon- 
tal length is about twice the critical wavelength 
(A c « 6.3z d ). The resolution is 32 x 32 x 32 so 
that one horizontal critical wavelength spans ~16 
cells, and one vertical scale length z d spans 4 cells. 
These choices for domain size and resolution are 
tested in §3.2.2. The simulation is terminated at 
~10.3f2 _1 , at which point the timestep has be- 
come three orders of magnitude smaller than the 
initial timestep (see the final paragraph of §2.1.1). 

Figure 5 displays a time series of the volume- 
rendered dust density in the bottom half of the 
box. Over the course of several dynamical times, 
density waves shear and amplify, eventually con- 
centrating into a single azimuthally elongated fila- 
ment. This filament then fragments radially. The 



fragments gravitationally scatter and merge; by 
the end of the simulation, two clumps remain. 

A simple diagnostic that we use throughout this 
paper is the time evolution of the maximum dust 
density, shown in the left panel of Figure 6. Com- 
parison with Figure 5 reveals that maxpd grows 
exponentially when the filament fragments radi- 
ally. The maximum dust density ceases to rise 
once the clumps finish coalescing. At this point 
each clump is gravitationally bound, with a max- 
imum central density that depends on the simula- 
tion resolution (§3.2.2). 
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Fig. 6. — Left: Time evolution of maximum dust 
density for run STD32 (= S). Right: Time evolution 
of kinetic energies averaged horizontally and vertically 
over a thin slab subtending two grid cells at the mid- 
plane (red = :r-component of kinetic energy; blue = y; 
green = z; black = total). 

The right panel of Figure 6 shows the time 
evolution of various kinetic energy densities, eval- 
uated in the three directions and excluding the 
background Keplerian shear. The energy densi- 
ties are averaged horizontally and vertically over 
a thin slab subtending two grid cells at the mid- 
plane (qualitatively similar results are obtained 
over larger vertical averages). The horizontal ki- 
netic energies grow exponentially from t — 2- 
7fi~ , with an exponential growth rate of ^1.5i7. 
Radial motions dominate azimuthal motions until 
the end of the simulation when they become com- 
parable. Vertical motions develop immediately af- 
ter the beginning of the simulation because our 
discretized initial conditions cannot be in perfect 
hydrostatic balance; however the magnitude of the 
vertical motions is small and stays roughly con- 
stant for t < 6f2 _1 . For t > 6fl~ 1 , vertical mo- 
tions amplify but for the most part remain smaller 
than horizontal motions. 

The in-plane motions of the dusty clumps are 
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Fig. 5. — Time evolution of gravitational instability in our standard 3D stratified dusty disk (run S = STD32). 
Shown are volume renderings of dust density for the bottom half of the disk at t = 0, 5.0, 6.4, 7.3, 8.0, and 10.3O -1 
(left to right, top to bottom). 



11 



Table 3 

3D Simulation Parameters ("Science Runs") 



Name 


Q s 


Mo 




(cg/n) 


2d 
(Cg/fi) 


(2d) 


Lx x 


Ly X 

(*S) 


L z 


Resolution 


Duration 

(n- 1 ) 


Qa 


Po 
(P f ) 


P& 
(P f ) 


Pi\ 
(P f ) 


Gl <c) 


S 


24 


35.0 


8.0 


0.20 


0.083 


6.3 


14 


< 14 > 


8 


32 > 


< 32 x 


32 


10.3 


0.5 


1.20 


0.16 


0.30 


Y 


Rl 


24 


165.0 


2.0 


0.91 


0.011 


42.6 


90 


< 90 > 


8 


256 ) 


< 256 


x 32 


11 


0.93 


1.21 


0.16 


1.05 


Y/N (<i) 


R2 


12 


93.0 


0.67 


1.04 


0.008 


150.3 


256 


< 256 


x 8 


256 > 


< 256 


x 32 


30 


1.86 


1.20 


0.16 


4.09 


N 


R3 


12 


143.0 


0.54 


1.07 


0.004 


255.4 


400 


< 400 


x 8 


400 > 


< 400 


x 32 


30 


1.86 


1.77 


0.16 


6.12 


N 


R4 


12 


322.0 


0.56 


1.07 


0.002 


220.7 


400 


< 400 


x 8 


400 > 


< 400 


x 32 


30 


1.2 


4.0 


0.16 


5.16 


N 


R5 


12 


322.0 


1.33 


0.87 


0.004 


44.7 


400 


< 400 


x 8 


400 > 


< 400 


x 32 


3.6 


0.5 


4.9 


0.16 


1.24 


Y 


SR 


24 


35.0 


8.0 


0.20 


0.083 


6.3 


400 


< 400 


x 8 


400 > 


< 400 


x 32 


10.0 


0.5 


1.20 


0.16 


0.30 


Y 


Z 


24 


35.0 


4.0 


0.52 


0.077 


13.1 


30 


< 30 > 


: 8 


32 > 


< 32 x 


32 


30 


1.0 


0.46 


0.16 


0.46 


N 


Q 


48 


35.0 


8.0 


0.35 


0.12 


8.6 


20 


< 20 > 


: 8 


32 > 


< 32 x 


32 


30 


1.0 


0.34 


0.16 


0.34 


N 


M 


24 


8.0 


8.0 


0.24 


0.71 


3.2 


8 


< 8 x 


8 


32 > 


< 32 x 


32 


30 


1.0 


0.25 


0.16 


0.25 


N 



Values are derived using Q* — 1. 

^ Values are derived using — 1. 

GI — Gravitational Instability. Y means max increases by orders of magnitude over a few dynamical times, and N means 
it does not. 

(d) Sec Figure 11. 
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Fig. 7. — Snapshot of the midplane for run S = 
STD32 at t = 9.6ft -1 . The largest in-plane velocity 
shown is 2.16 c g . 

illustrated in Figure 7 with a snapshot of the mid- 
plane slice of STD32 at t = 9.6 fi -1 . The dust 
clumps are seen spinning about their centers of 
mass as a consequence of angular momentum con- 
servation. 

3.2.2. Resolution and box size 

Table 4 lists the parameters of experiments de- 
signed to test our choices for resolution, box size, 



and grid-cell aspect ratio. 

Figure 8 shows how varying the resolution 
changes the evolution of our standard, gravita- 
tionally unstable run (STD32 — also labeled S 
in Table 3). We use again the simple metric of 
max vs. t. Broadly speaking, the runs STD16, 
STD32, STD64 are all "acceptable" insofar as they 
all yield increases in max p^ by orders of magni- 
tude within several dynamical times (t < 8f2 _1 ). 
By contrast, the lowest resolution run, STD8, is 
unacceptable. Thus, the minimum acceptable res- 
olution appears to be ^2 cells per scale length z& in 
the vertical direction (cf. Nelson 2006 who found 
that a minimum of four smoothing lengths per 
scale height is required for SPH simulations), and 
^8 cells per critical wavelength A c in the horizon- 
tal directions. Our standard choices for resolution 
— as well as the resolutions characterizing all our 
"science" runs, listed in Table 3 and discussed in 
§3.2.3 — satisfy these minimum requirements by 
a safety factor of 2. 

Examining Figure 8 more critically, we see that 
the maximum value attained by max has not 
converged with resolution. Increasing the resolu- 
tion enables us to resolve ever higher densities in 
the collapsing clumps. Another point of concern 
is the non-uniform aspect ratios of individual grid 
cells, which ranges from x:y:z w 2:2:1 to 4:4:1 over 
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Table 4 

3D Simulation Parameters to Test Box Size and Resolution 



Name 


L x x Ly x L z (z^ L ) 


Resolution 




Duration (Q ) 


STD32(S) 


14 x 14 x 8 


32 x 32 x 32 


Y 


9.8 


STD8 


14 x 14 x 8 


8x8x8 


N 


30.0 


STD16 


14 x 14 x 8 


16 x 16 x 16 


Y 


11.0 


STD64 


14 x 14 x 8 


64 x 64 x 64 


Y 


11.0 


U32 


14 x 14 x 8 


56 x 56 x 32 


Y 


10.0 


LZ2 


14 x 14 x 2 


32 x 32 x 8 


N 


30.0 


LZ4 


14 x 14 x 4 


32 x 32 x 16 


Y 


8.5 


LZ6 


14 x 14 x 6 


32 x 32 x 24 


Y 


8.0 


LZ10 


14 x 14 x 10 


32 x 32 x 40 


Y 


9.0 


LZ14 


14 x 14 x 14 


32 x 32 x 56 


Y 


10.5 


LXY6 


6x6x8 


16 x 16 x 32 


N 


30.0 


LXY10 


10 x 10 x 8 


24 X 24 X 32 


Y 


8.6 


LXY20 


20 x 20 x 8 


48 x 48 x 32 


Y 


8.7 



' a ' GI = Gravitational Instability. Y means max p<j increases by orders of 
magnitude over a few dynamical times, and N means it does not. 
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Fig. 8. — Time evolution of the maximum dust den- 
sity for the resolution study (see Table 4). 



Fig. 9. — Time evolution of the maximum dust den- 
sity for our box size tests (see Table 4). 



our set of science simulations (Table 3). The run 
U32 is characterized by perfectly cubical grid cells 
(1:1:1); the evolution is similar to STD32, but is 
characterized by an earlier onset of gravitational 
instability, and stronger density fluctuations. This 
comparison suggests that our science runs with 
non-cubical grid cells are biased slightly against 
gravitational instability. 

We next investigate how box size affects our re- 
sults. For all box size experiments, the spatial res- 
olution is kept at its standard value (32 grid cells 
per 14zd in either horizontal direction, and 4 grid 



cells per Zd in the vertical direction). Runs LZ2 
through LZ14 vary box height L z while keeping 
L x and L y fixed at their standard (STD32 = S) 
values. As Figure 9 reveals, box heights of 4-14^ 
yield comparable results, while a box height of 2z^ 
is unacceptable. For the most part, increasing the 
box height seems to delay the onset of gravita- 
tional instability, with LZ4 being the exception to 
this rule. 

Our 2D simulations indicated that L x and L y 
must be large enough to encompass at least one 
critical wavelength A c . Our 3D simulations bear 
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out this same requirement. Figure 9 shows that 
run LXY6, for which the box size is just under one 
critical wavelength, does not exhibit gravitational 
instability, unlike its bigger box counterparts. 

To summarize our findings in this subsection: 
(1) The simulation box should be at least 4zd 
tall (2zd above and below the midplane). (2) 
Each horizontal dimension must be longer than 
one critical wavelength A c as given by equation 
(30). (3) Simulations require a vertical resolu- 
tion of > 2 grid cells per scale length Zd, and a 
horizontal resolution of > 8 grid cells per critical 
wavelength. (4) Individual grid cells that have in- 
creasingly non-uniform aspect ratios (squatter ver- 
tically than horizontally) tend to suppress gravita- 
tional instability, but the bias is minor and aspect 
ratios up to 4:4:1 appear acceptable. All of our 
science simulations (Table 3; §3.2.3) satisfy these 
requirements, in some cases by factors of 2. 

3.2.3. Criteria for gravitational collapse 

Table 3 lists the simulations designed to test 
which of the various proposed criteria for gravita- 
tional instability is the best predictor of collapse. 
Figures 2 and 10 describe the initial dust and gas 
profiles, while Figure 11 displays the results using 
our simple diagnostic of max pd vs. time. 

First consider runs S and R1-R5, and ask 
whether these runs favor p\ or pjj for the density 
required for gravitational collapse. Because dust 
is a major component of our disks, we do not ex- 
pect p\ — which is strictly valid only for pure gas 
disks — to be a good predictor. Indeed in all six of 
these runs, the midplane density po exceeds p*, by 
factors of 7.5-30, yet only runs S and R5, and to 
a much lesser extent Rl, exhibit collapse. All six 
runs indicate instead that pjj — equivalently, 
— is the better predictor, with the critical value 

0.5<Q d <0.9. (33) 

There is some concern that the comparison be- 
tween runs R2-R5 and run S may not be fair be- 
cause runs R2-R5 have a factor of ~2 poorer spa- 
tial resolution in x and y compared to run S. This 
concern is allayed by run SR, which has the same 
physical parameters as S but is run with the box 
size and resolution of R3, and which turns out to 
behave qualitatively similarly to S (see Figure 11). 

Our conclusion that Pjj is relevant and that Q* d 
obeys (33) is supported further by runs Z, Q, and 
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Fig. 10. — Initial conditions for our science runs 
which explore parameter space. Solid lines denote 
dust, and dashed lines denote gas. The vertical lines 
delimit the vertical boundaries of our simulation box. 

M, each of which varies one of the three input 
parameters Sd/S g , Qd, and po. 

Although runs R2-R4 do not exhibit the dra- 
matic growth in p& shown by runs S, SR, and R5 
— a result that we interpret to mean that p"^ gives 
the correct criterion for gravitational collapse - 
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Fig. 11. — Time evolution of the maximum dust den- 
sity in our science simulations. Only for runs S, SR, 
and R5 does po > pfi, and indeed only those runs 
exhibit dramatic growth of the dust density due to 
gravitational instability. 




Fig. 12. — Snapshots of the dust density at the mid- 
plane for run R3 (left panel) and run SR (right panel) . 
In R3, the initial midplane density po > Psckiyai an d 
there is some clumping, but it is much lower in ampli- 
tude compared to run SR, for which po > pj r . 

runs R2-R4 do show some clumping. Figure 12 
compares snapshots of runs R3 and SR (performed 
with the same box size and resolution), taken at 
the same time t — 10f2 _1 . Filaments do form in 
R3, although they are much weaker in density con- 
trast compared to the filaments in SR. The mild 
growth shown in runs R2 and R3 might simply re- 
flect the fact that their values for Qd = 1.86 are 
still too close to Q d to suppress instability entirely. 
An alternative (and not mutually exclusive) pos- 
sibility is that because p > Pg ki ya = 0.60p^ for 
runs R2-R4, the disk might be exhibiting the un- 
stable (and formally incompressible) mode found 
by Sekiya (1983). Whatever the interpretation, 
the modest growth factors exhibited by R2-R4 



seem unlikely to lead to planetesimal formation. 

Finally, what about p^ ochc vs. pjj? Here runs 
R4 and R5 are the most telling. Both runs are 
characterized by the largest midplane densities 

Po > PRocho but onl y R5 > for which po > 
undergoes gravitational collapse (see Figure 11). 

Table 5 summarizes how the various candi- 
date critical densities relate to one another and to 
the midplane density for our science simulations. 
From Table 5, pjj emerges as the best predictor of 
collapse. 

4. SUMMARY AND DISCUSSION 




Fig. 13. — A tale of two particle sublayers, one of 
which is thinner and denser than the other. Dust den- 
sity is plotted as a solid line, and gas density as a 
dashed line. The disks have identical masses and bulk 
metallicities, enhanced over those of the minimum- 
mass solar nebula by factors of 3-4. Left: Midplane 
density p — Pr oc1ic = 3.5pt and Qd = 10.4. Right: 
Midplane density po « Pn w 10 2 /0r oc1io and Qd = 1. 
According to the results of our simulations, only the 
model in the right panel, having the thinner and denser 
sublayer, should be on the verge of gravitational col- 
lapse — in the limit that particles are aerodynamically 
perfectly coupled to gas. We argue in §4.1 that when 
the perfect coupling approximation breaks down, it 
may be possible for the disk on the left to undergo 
gravitational instability. 

Dust grains settle toward the midplanes of pro- 
toplanetary disks, forming a sublayer of solid par- 
ticles sandwiched from above and below by gas. 
Whether this sublayer can become thin enough 
and dense enough to undergo gravitational insta- 
bility and fragment into planetesimals is an out- 
standing question. We have found in this work 
that the density threshold for gravitational col- 
lapse can be extraordinarily high — much higher 
even than the Roche density PR 0che = 3.5M*/?- 3 , 



15 



Table 5 

Comparison of Critical Densities and Actual Midplane Density for Science Simulations 



Name 


Critical density relations 






S 


Pi < Ph. < P&kiya < PO < PRochc 




Y 


HI 


Pi < P&kiya < PlI < PO < P R ochc 




Y/N( 6 ) 


R2 


Pi < PSckiya < PO < P Roc hc 


< p r u 


N 


R3 


Pi < PSckiya < PO < PRochc 


< ph 


N 


R4 


Pi* < P S ckiya < PRochc < P° < P*U 




N 


R5 


Pi < PSckiya < Ph < PRochc < PO 




Y 



( a ) GI = Gravitational Instability. Y means maxpj increases by 
orders of magnitude over a few dynamical times, and N means it does 
not. 

W See Figure 11. 



where M* is the mass of the central star and r is 
the orbital radius. To trigger collapse in the limit 
that dust particles are small enough to be tightly- 
coupled to gas, the density po in the sublayer must 
be such that the Toomre stability parameter 



Pu 
Po 



1/2 



< 1 



(34) 



(35) 



where 

(For more precise relations, see equations 8, 13, 
and 33.) Here Q g is the Toomre parameter for 
the ambient (and much thicker) gas disk, Ed/S g 
is the ratio of surface densities of dust and gas 
(i.e., the height-integrated metallicity), and 0.5 < 
Q d < 0.9 as measured from our simulations. For 
an astrophysically plausible disk having 3x the 
mass of the minimum-mass solar nebula (Q g « 10) 
and a bulk metallicity enriched over solar by a 
factor of 4 (Sd/S g s» 0.06), the critical density 



p* u « 1.3 x 10 2 Qr 2 PRochc 



(36) 



Figure 13 portrays two sublayers — one for 
which po = p^ ochc and another, much thinner sub- 
layer for which p « p* u « 10 2 p^ ocho (Q d = 1). 
The results of our simulations, performed in the 
limit of perfect aerodynamic coupling between 
particles and gas, indicate that only the latter, 
much denser disk is on the verge of fragmenting. 

Qualitatively, such extraordinary densities are 
required for gravitational instability because gas 



pressure renders the sublayer extremely stiff. 
Sound-crossing times for thin layers are easily 
shorter than free-fall times. We can examine the 
competition between stabilizing pressure, stabi- 
lizing rotation, and de-stabilizing self-gravity in 
both the horizontal (in-plane) and vertical direc- 
tions. Horizontal stability is controlled by Qd- 
when Qd > Q d ~ 1, all horizontal lengthscales 
A < 2cj/GSd are stabilized by pressure, and all 
scales A > 2c 2 /GSd are stabilized by rotation, 
where Cd is the effective sound speed in the dust- 
gas mixture. At the same time, vertical stabil- 
ity is assured whenever the sound-crossing time 
across the vertical thickness of the sublayer 2Hd 
is shorter than the free-fall time: 



2F d I 



VG Pd 



(37) 



which, after substituting Hd ~ Sd/2/Od and Cd 
c g \fp~gfpA, translates to 



1 7T 



(38) 



which is easily satisfied for reasonable disk param- 
eters. 

Taken at face value, the higher density thresh- 
old pjj that our work establishes argues against 
using aerodynamically well-coupled particles to 
form planetesimals. The Kelvin-Helmholtz insta- 
bility probably stops the midplane density from 
approaching p^. For disks enriched in both mass 
and metallicity (by factors of several over the 
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minimum- mass solar nebula), midplane densities 
might just barely reach ~M*/r 3 , but probably 
will not exceed it (e.g., Chiang 2008, Lee et al. 
2010a,b). 2 

4.1. Limitations of Our Study, and 
Directions for Future Research 

We have worked in the limit that the stopping 
times £ s top of particles in gas are small compared to 
all other timescales. But in reality, finite particle 
sizes imply finite t s t op (see Figure 1). When the 
assumption of infinitesimal stopping time breaks 
down, new effects may appear that might lower 
the threshold for gravitational instability. 

One such effect is as follows. Consider again 
the competition between stabilizing pressure and 
de-stabilizing self-gravity (in either the vertical or 
horizontal directions). A major reason why the 
sublayer so strongly resists collapse is that sound 
waves travel quickly across it. We have taken 
the sound speed for our dust-gas suspension to 
be Cd = c g / ^/l + pd/Pg ~ C g\/Pg/Pd (equations 9 
and 10). But this presumes that particles are per- 
fectly coupled to gas. If the sound-crossing time 
across some scale A were to become shorter than 
the particle stopping time, i.e., if 

A A 



< t 



stop 



(39) 



then our use of Cd ~ Cg-y/Pg/Pd would be invalid. 
Particles on scales A would lose support from gas 
pressure and become susceptible to gravitational 
instability. 

To get a sense of where in parameter space this 
instability may lie, we normalize A to the full ver- 
tical thickness of the sublayer: 



A = 2iJ ri A 



ASd 
Pd 



(40) 



where A can take any value (larger than or smaller 
than unity). Then equation (39) for the loss of 
pressure support translates to a midplane density 
(dominated by dust) of 

2 



PO ~ Pd 



> 



2 M* 



7r r° 



X 2 



1 



s g ; Qg(nt stop y 



(41) 



2 See, however, Sekiya (1998), who found that self- 
gravitating, non-rotating sublayers having constant 
Richardson number could develop cusps of infinite density 
at the midplane. 



where fl is the Kepler orbital frequency. For self- 
gravity to resist tidal disruption, pd = PRochc = 
3.5M*/r 3 . Substituting this requirement into 
(41), we find that 



^^stop ^ 




for particles on scales A to decouple from sound 
waves. For A = 1, requirement (42) could be ful- 
filled by particles having sizes of a few millimeters 
to a few centimeters at distances of 1-10 AU (Fig- 
ure 1 — but note that the curves in the figure need 
to be adjusted by factors of a few for mass-enriched 
nebulae). For A < 1, even smaller particles could 
lose pressure support and collapse gravitationally. 

Future simulations that include finite particle 
stopping times could try to find such an instabil- 
ity. A complication would be that accounting for 
finite t st0 p would introduce the streaming insta- 
bility, which could prevent the dust density from 
attaining the Roche value — see, e.g., runs R21-3D 
and R41-3D in Figure 5 of Bai & Stone (2010), for 
which fttstop < 0.1 and pd < /Or oc i 1c - To find the 
instability that we are envisioning, one would have 
to restrict ftt s to P to small enough values to sup- 
press the streaming instability — thereby permit- 
ting the setting of grains into sublayers for which 
Pd = /°Rochc — while at the same time keeping 
fiistop large enough to satisfy (42) and nullify pres- 
sure support. 

We are grateful to Eve Ostriker and Andrew 
Youdin for discussions. Section 4.1 was inspired by 
discussions with Eve that clarified the limitations 
of our study and pointed the way to a possible 
new route to gravitational instability. We thank 
Xucning Bai, Chang-Goo Kim, Eve Ostriker, Ian 
Parrish, and Jim Stone for help in augmenting 
Athena. The simulations were performed with the 
Berkeley cluster Henyey, which was made possi- 
ble by a National Science Foundation Major Re- 
search Instrumentation (NSF MRI) grant. Finan- 
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